rm(list = ls())

wd <- ".../Replication/"
setwd(wd)

# Load/install packages --
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  foreign, 
  ggplot2, 
  estimatr,
  texreg,
  xtable,
  fastDummies,
  sandwich,
  dplyr,
  janitor, 
  gridExtra,
  gsheet,
  zoo,
  interflex,
  lubridate,
  tidyverse,
  stringi,
  readxl,
  ri2,
  modelsummary,
  ggpubr
)

options(scipen=999)

# Note: given the random component inherent to the randomization tests performed
# by ri(), you may observe small differences in the p-values of the randomization
# tests performed using that function.

# Load data ----
data <- read.csv("Data/baseline database.csv")
data$imprisoned_harvest <- (data$dec+data$jan+data$feb)>0
data$imprisoned_noharvest <- (data$mar+data$apr+data$may+data$jun+data$jul+data$aug+data$sep+data$oct+data$nov)>0

censo <- read.csv("Data/census_nonwhites.csv")

# select data from winners
winners <- subset(data, data$group=="Liberado sorteo")
win <- winners %>% select(imprisoned_dummy, 
                          lottery_winner,
                          imprisoned_harvest,
                          imprisoned_noharvest,
                          robbery, 
                          assault, 
                          discretionary)

names(win)[names(win)=="full_name"] <- "name"

# clean data from census
censo$lottery_winner <- 0
censo$imprisoned_harvest <- as.numeric(grepl("december", censo$months) | grepl("january", censo$months) | grepl("february", censo$months))
censo$imprisoned_noharvest <- as.numeric(grepl("march", censo$months) | grepl("april", censo$months) | grepl("may", censo$months) | 
                                                           grepl("june", censo$months) | grepl("july", censo$months) | grepl("august", censo$months) |
                                                           grepl("september", censo$months) | grepl("october", censo$months) | grepl("november", censo$months))

# remove from the census the individuals who might be the winners
censo$name1 = paste(censo$nombre,censo$apellido)
censo$name2 = paste(censo$nombre,censo$apellido_jefehogar)
censo$name_search = ifelse(is.na(censo$apellido)==TRUE,censo$name2,censo$name1)
censo$name_search = ifelse(is.na(censo$apellido)==TRUE & is.na(censo$apellido_jefehogar)==TRUE,NA,censo$name_search)
censo = subset(censo, is.na(censo$name_search)==FALSE)


possiblematches = NA

for(i in 1:nrow(winners)){
  list = agrep(winners$full_name[i], 
                censo$name_search, 
                max.distance =0.10, value = TRUE)
  possiblematches = c(possiblematches,list)
  print(i)
}

censo2 <- subset(censo, !(censo$name_search %in% possiblematches))

# create different subsets of enslaved persons in the census
enslaved <- subset(censo2, censo2$slave==1)
names(enslaved)[names(enslaved)=="imprisoned"] <- "imprisoned_dummy"

# create subsets of the census and list of winners

all_15_20 <- enslaved %>% filter(age>=15 & age<=20) %>% 
  select(imprisoned_dummy, lottery_winner,
         imprisoned_harvest,
         imprisoned_noharvest,
         robbery, 
         assault, 
         discretionary)

all_15_30 <- enslaved %>% filter(age>=15 & age<=30) %>% 
  select(imprisoned_dummy, lottery_winner,
         imprisoned_harvest,
         imprisoned_noharvest,
         robbery, 
         assault, 
         discretionary)

all_15_40 <- enslaved %>% filter(age>=15 & age<=40) %>% 
  select(imprisoned_dummy, lottery_winner,
         imprisoned_harvest,
         imprisoned_noharvest,
         robbery, 
         assault, 
         discretionary)

all <- enslaved %>% select(imprisoned_dummy, lottery_winner,
                           imprisoned_harvest,
                           imprisoned_noharvest,
                           robbery, 
                           assault, 
                           discretionary)


# final data for analyses
d1 <- rbind.data.frame(win,all)
d2 <- rbind.data.frame(win,all_15_20)
d3 <- rbind.data.frame(win,all_15_30)
d4 <- rbind.data.frame(win,all_15_40)


# compute exact p-values
declare.d1 <- declare_ra(N = nrow(d1), m=sum(d1$lottery_winner)) 
declare.d2 <- declare_ra(N = nrow(d2), m=sum(d2$lottery_winner)) 
declare.d3 <- declare_ra(N = nrow(d3), m=sum(d3$lottery_winner)) 
declare.d4 <- declare_ra(N = nrow(d4), m=sum(d4$lottery_winner)) 




#### TABLE A2: Effect of emancipation, harvest vs. nonharvest months

# Harvest months
all_harvest<-summary(conduct_ri(formula = imprisoned_harvest ~ lottery_winner,
                                declaration = declare.d1, sharp_hypothesis = 0, assignment = "lottery_winner",
                                data = d1, sims = 10000))
all_15_20_harvest<-summary(conduct_ri(formula = imprisoned_harvest ~ lottery_winner,
                                      declaration = declare.d2, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d2, sims = 10000))
all_15_30_harvest<-summary(conduct_ri(formula = imprisoned_harvest ~ lottery_winner,
                                      declaration = declare.d3, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d3, sims = 10000))
all_15_40_harvest<-summary(conduct_ri(formula = imprisoned_harvest ~ lottery_winner,
                                      declaration = declare.d4, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d4, sims = 10000))

# Non-harvest months
all_non_harvest<-summary(conduct_ri(formula = imprisoned_noharvest ~ lottery_winner,
                                    declaration = declare.d1, sharp_hypothesis = 0, assignment = "lottery_winner",
                                    data = d1, sims = 10000))
all_15_20_non_harvest<-summary(conduct_ri(formula = imprisoned_noharvest ~ lottery_winner,
                                          declaration = declare.d2, sharp_hypothesis = 0, assignment = "lottery_winner",
                                          data = d2, sims = 10000))
all_15_30_non_harvest<-summary(conduct_ri(formula = imprisoned_noharvest ~ lottery_winner,
                                          declaration = declare.d3, sharp_hypothesis = 0, assignment = "lottery_winner",
                                          data = d3, sims = 10000))
all_15_40_non_harvest<-summary(conduct_ri(formula = imprisoned_noharvest ~ lottery_winner,
                                          declaration = declare.d4, sharp_hypothesis = 0, assignment = "lottery_winner",
                                          data = d4, sims = 10000))


Table_A2<-data.frame(matrix(nrow = 8, ncol = 4))
colnames(Table_A2)<-c('Control Group','','Harvest Months','Nonharvest Months')
Table_A2[1,1]<-'15 and older'
Table_A2[2,1]<-'15 and older'
Table_A2[3,1]<-'15-40'
Table_A2[4,1]<-'15-40'
Table_A2[5,1]<-'15-30'
Table_A2[6,1]<-'15-30'
Table_A2[7,1]<-'15-20'
Table_A2[8,1]<-'15-20'

Table_A2[c(1,3,5,7),2]<-'Diff. in Means'
Table_A2[c(2,4,6,8),2]<-'ri p-value'

Table_A2[1,3]<-all_harvest[,'estimate']
Table_A2[2,3]<-all_harvest[,'two_tailed_p_value']
Table_A2[3,3]<-all_15_40_harvest[,'estimate']
Table_A2[4,3]<-all_15_40_harvest[,'two_tailed_p_value']
Table_A2[5,3]<-all_15_30_harvest[,'estimate']
Table_A2[6,3]<-all_15_30_harvest[,'two_tailed_p_value']
Table_A2[7,3]<-all_15_20_harvest[,'estimate']
Table_A2[8,3]<-all_15_20_harvest[,'two_tailed_p_value']

Table_A2[1,4]<-all_non_harvest[,'estimate']
Table_A2[2,4]<-all_non_harvest[,'two_tailed_p_value']
Table_A2[3,4]<-all_15_40_non_harvest[,'estimate']
Table_A2[4,4]<-all_15_40_non_harvest[,'two_tailed_p_value']
Table_A2[5,4]<-all_15_30_non_harvest[,'estimate']
Table_A2[6,4]<-all_15_30_non_harvest[,'two_tailed_p_value']
Table_A2[7,4]<-all_15_20_non_harvest[,'estimate']
Table_A2[8,4]<-all_15_20_non_harvest[,'two_tailed_p_value']

View(Table_A2)

print(xtable(Table_A2), file = "Output/Table A2.tex")

#### TABLE A3: Effect of emancipation by type of criminal charge

# Robbery
all_robbery<-summary(conduct_ri(formula = robbery ~ lottery_winner,
                                declaration = declare.d1, sharp_hypothesis = 0, assignment = "lottery_winner",
                                data = d1, sims = 10000))
all_15_20_robbery<-summary(conduct_ri(formula = robbery ~ lottery_winner,
                                      declaration = declare.d2, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d2, sims = 10000))
all_15_30_robbery<-summary(conduct_ri(formula = robbery ~ lottery_winner,
                                      declaration = declare.d3, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d3, sims = 10000))
all_15_40_robbery<-summary(conduct_ri(formula = robbery ~ lottery_winner,
                                      declaration = declare.d4, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d4, sims = 10000))

# Assault
all_assault<-summary(conduct_ri(formula = assault ~ lottery_winner,
                                declaration = declare.d1, sharp_hypothesis = 0, assignment = "lottery_winner",
                                data = d1, sims = 10000))
all_15_20_assault<-summary(conduct_ri(formula = assault ~ lottery_winner,
                                      declaration = declare.d2, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d2, sims = 10000))
all_15_30_assault<-summary(conduct_ri(formula = assault ~ lottery_winner,
                                      declaration = declare.d3, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d3, sims = 10000))
all_15_40_assault<-summary(conduct_ri(formula = assault ~ lottery_winner,
                                      declaration = declare.d4, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d4, sims = 10000))

# Discretionary
all_discretionary<-summary(conduct_ri(formula = discretionary ~ lottery_winner,
                                      declaration = declare.d1, sharp_hypothesis = 0, assignment = "lottery_winner",
                                      data = d1, sims = 10000))
all_15_20_discretionary<-summary(conduct_ri(formula = discretionary ~ lottery_winner,
                                            declaration = declare.d2, sharp_hypothesis = 0, assignment = "lottery_winner",
                                            data = d2, sims = 10000))
all_15_30_discretionary<-summary(conduct_ri(formula = discretionary ~ lottery_winner,
                                            declaration = declare.d3, sharp_hypothesis = 0, assignment = "lottery_winner",
                                            data = d3, sims = 10000))
all_15_40_discretionary<-summary(conduct_ri(formula = discretionary ~ lottery_winner,
                                            declaration = declare.d4, sharp_hypothesis = 0, assignment = "lottery_winner",
                                            data = d4, sims = 10000))


Table_A3<-data.frame(matrix(nrow = 8, ncol = 5))
colnames(Table_A3)<-c('Control Group','','Robbery','Assault','Discretionary')
Table_A3[1,1]<-'15 and older'
Table_A3[2,1]<-'15 and older'
Table_A3[3,1]<-'15-40'
Table_A3[4,1]<-'15-40'
Table_A3[5,1]<-'15-30'
Table_A3[6,1]<-'15-30'
Table_A3[7,1]<-'15-20'
Table_A3[8,1]<-'15-20'

Table_A3[c(1,3,5,7),2]<-'Diff. in Means'
Table_A3[c(2,4,6,8),2]<-'ri p-value'

Table_A3[1,3]<-all_robbery[,'estimate']
Table_A3[2,3]<-all_robbery[,'two_tailed_p_value']
Table_A3[3,3]<-all_15_40_robbery[,'estimate']
Table_A3[4,3]<-all_15_40_robbery[,'two_tailed_p_value']
Table_A3[5,3]<-all_15_30_robbery[,'estimate']
Table_A3[6,3]<-all_15_30_robbery[,'two_tailed_p_value']
Table_A3[7,3]<-all_15_20_robbery[,'estimate']
Table_A3[8,3]<-all_15_20_robbery[,'two_tailed_p_value']

Table_A3[1,4]<-all_assault[,'estimate']
Table_A3[2,4]<-all_assault[,'two_tailed_p_value']
Table_A3[3,4]<-all_15_40_assault[,'estimate']
Table_A3[4,4]<-all_15_40_assault[,'two_tailed_p_value']
Table_A3[5,4]<-all_15_30_assault[,'estimate']
Table_A3[6,4]<-all_15_30_assault[,'two_tailed_p_value']
Table_A3[7,4]<-all_15_20_assault[,'estimate']
Table_A3[8,4]<-all_15_20_assault[,'two_tailed_p_value']

Table_A3[1,5]<-all_discretionary[,'estimate']
Table_A3[2,5]<-all_discretionary[,'two_tailed_p_value']
Table_A3[3,5]<-all_15_40_discretionary[,'estimate']
Table_A3[4,5]<-all_15_40_discretionary[,'two_tailed_p_value']
Table_A3[5,5]<-all_15_30_discretionary[,'estimate']
Table_A3[6,5]<-all_15_30_discretionary[,'two_tailed_p_value']
Table_A3[7,5]<-all_15_20_discretionary[,'estimate']
Table_A3[8,5]<-all_15_20_discretionary[,'two_tailed_p_value']

View(Table_A3)

print(xtable(Table_A3), file = "Output/Table A3.tex")


